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ABSTRACT 

We derive analytical expressions for the change in the orbital eccentricity of a binary 
following a distant encounter with a third star on a hyperbolic or parabolic orbit. 
To establish the accuracy of these expressions, we present detailed comparisons with 
the results of direct numerical integrations of the equations of motion for the three 
bodies. We treat with particular care the difficult case of a binary with zero initial 
eccentricity. In this case, we show that the eccentricity Se induced by the encounter 
declines in general as a power-law, 8e oc (a/r p ) 5 / 2 , where a is the binary semi-major 
axis and r p is the periastron distance of the encounter. This power-law arises from the 
octupole- level secular perturbation of the binary. In contrast, non-secular quadrupole- 
level perturbations induce an eccentricity change that declines exponentially with r p . 
These non-secular effects can become dominant at sufficiently small r p , for a suffi- 
ciently high relative velocity, or for a sufficiently massive perturber. We also derive 
cross sections for eccentricity change and compare our results with those of previous 
studies based on numerical scattering experiments. Our results have important impli- 
cations for a number of astrophysical problems including, in particular, the evolution 
of binary millisecond pulsars in globular clusters. 
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1 INTRODUCTION 

Three-body scattering plays an essential role in stellar dynamics. In particular the interactions between binary stars and 
single stars in globular clusters are thought to be an important dynamical process, providing support against core collapse 
(Spitzer 1987; Hut et al. 1992). What is important here is the energy exchanged, but there are applications in which it is 
also important to understand the changes in the eccentricity of a binary caused by an encounter with another star. In some 
cases even very small changes in eccentricity can be very significant. Most importantly, binary millisecond pulsars have been 
discovered in large numbers in globular clusters and their orbital eccentricities can be measured with extraordinary precision. 
Although the orbits of binary millisecond pulsars in the galactic plane are nearly circular (for example, PSR J2317+1439 
has a measured eccentricity of 1.2 x 10 -6 ; see Camilo et al. 1993), those found in the dense stellar environment of globular 
clusters may be perturbed by the high rate of encounters with other stars. If, therefore, the rate at which the eccentricity is 
changed by encounters is known, interesting constraints can sometimes be placed on the lifetimes of the binaries (Rappaport 
et al. 1989; Rasio & Heggie 1995). 

Much of our knowledge of the effect of encounters on eccentricity is numerical (Hills 1975; Hills & Fullerton 1980; Hut 
& Paczyhski 1984; Rappaport et al. 1989). In the present paper the approach is analytical, with numerical work used only 
to confirm the theoretical results. Heggie (1975) gave a theoretical estimate (based on first-order perturbation theory) for the 
change in eccentricity, and found that it varies inversely as the 3/2 power of r p , the distance of closest approach of the third 
body. Unfortunately, his result vanishes in the important case in which the initial orbit is circular. Phinney (1992) and Heggie 
(unpublished) then found, again analytically, that in this case the induced eccentricity appeared to decrease exponentially 
with a power of r p , though the results were restricted to special situations (e.g., equal masses, coplanar motions, and parabolic 
or rectilinear motion of the third body). While the authors of the present paper were reworking the theory for more general 
cases, S.L.W. McMillan pointed out to us, on the basis of numerical experiments with a new scattering package (McMillan 
1996), that the predicted exponential dependence for initially circular binaries fails when the masses of the binary components 
are unequal (which clearly is the important case in applications). This discovery stimulated much of the remaining work for 



2 D. C. Heggie and F. A. Rasio 



the present paper. 

In this paper our aim is to clarify the regimes in which the various theoretical results are applicable, and to extend them, 
where necessary, to arbitrary masses and arbitrary geometries. The one common assumption is that the encounter is tidal, 
i.e., that r p o, where a is the initial semi-major axis of the binary, though it may be necessary to take account of the masses 
of the participants in this condition when the masses are very different. We also usually assume that the encounter is slow. 

In the following section we outline the theoretical results, but leave detailed derivations to an appendix. We also compare 
the theoretical results with new numerical data, designed to explore the full dependence on r p , the masses, the initial eccen- 
tricity of the binary, and the geometry of the orbits. In Section 3 we convert our results into cross sections, which are then 
compared with results of other authors. In the last section we summarise our conclusions and outline briefly some applications. 
Another summary, which places particular emphasis on applications to binary pulsars in globular clusters, can be found in 
Rasio & Heggie (1995). 



2 GENERAL EXPRESSIONS FOR THE CHANGE IN ECCENTRICITY 
2.1 Notations, Assumptions and Basic Equations 

After setting out our notation and assumptions, in this section we shall give the basic equations for the perturbation of a 
binary by a third body. Much of this is elementary, but it is helpful because it makes some properties of our detailed results 
quite obvious. 

We consider an encounter between a single star of mass mz and a binary whose components have masses mi and m-2, 
and we define M12 = mi + m,2 and M123 = M12 + m^. Long before the encounter begins, the eccentricity of the binary is e, 
(though often we shall assume this is zero), and its semi-major axis is a. The orbit of the third body relative to the barycentre 
of the binary is hyperbolic with eccentricity e'. We denote its impact parameter by p and its speed relative to the barycentre, 
while they are still far apart, by V. Then we have 



p 2 V A 

and we shall be able to obtain results for a parabolic encounter by letting V — » 0, or e' — » 1. 

We make two assumptions about the encounter, that it is both tidal and slow. For the first of these, we assume that the 
closest distance of approach of the third body, which will be denoted by r p , is still considerably larger than a. The ratio of 
these distances will play the role of a perturbation parameter. For the second assumption, we require that the angular velocity 
of the perturber, when it is closest to the binary, is considerably smaller than the angular velocity of the binary components. 
It is easily shown that this is equivalent to the condition 



M12 /r p \ 3 / 2 



, (M'" * »1, (2) 
M123 V a ) v^TT 

and so the second assumption is implied by the first unless the masses are very disparate or e' is very large. (This last case, 
the "impulsive limit", is briefly considered in Section 3.2 and Appendix A4.) 

Now we develop the equations of motion. Let r be the position of m 2 relative to mi, and R the position of m 3 relative 
to the barycentre of the binary. Then the equation for the relative motion of the binary components is 

»-^ + ». (3) 

with 



Gm 3 - (-m 2 ) 

¥ -—2^ Vr 



(4) 



where P n is the nth Legendre polynomial. 

Except for a factor — mim 2 /Mi2, equation (4) is simply the gradient of the mutual potential of the components of the 
binary and the third body, expressed as a multipole expansion in powers of r/R. The monopole and dipole terms (i.e., those 
given by n = and 1) make no contribution to F, and will play no further part in this paper. 

Our aim is to compute the change in eccentricity of the binary over the entire duration of the encounter with the third 
body. To do this we have found it convenient to use the "eccentric axis" (Pollard 1976), which is a vector also associated with 
the names of Laplace, Hamilton, Lenz and Runge (Goldstein 1980). For unperturbed Keplerian motion it is a constant vector 
which points along the major axis of the elliptic orbit of the binary components, and its magnitude is e. Its definition is 

1 r 

-r x (r x r) . (5) 
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Using equation (3) it is easy to show that, in the presence of a perturbing acceleration F, its variation is determined by 

GM 12 e = 2(F ■ r)r — (r ■ r)F — (F ■ r)r (6) 

(cf. also Burdet 1967). These equations are the basis of the theory described in Appendix A, but some of the results we shall 
now summarise are more easily understood with reference to these formulae. 



2.2 Binaries with Non-Zero Initial Eccentricity 

We begin with what would seem to be the most straightforward and general case, that of a binary of arbitrary initial eccentricity 
e. Also, since we are assuming that r p >• a, it is natural to neglect all but the lowest non-vanishing term in equation (4). 
From what has been said, this is the quadrupole term n = 2 in the interaction between the binary and the third body. 

From the assumption that the encounter is slow (§2.1) it follows that, to lowest order, perturbations may be computed 
by first averaging over the fast relative motion of the components of the binary. The resulting estimate of the change in 
eccentricity is given (with a sign error) in the case of parabolic motion of the third body by Heggie (1975, eq. [5.66]), and it 
is rederived in the general case (and with the correct sign) in Appendix A to the present paper. 

In order to state this result, we introduce standard angles describing the relative orientation of the orbits. The line of 
nodes is defined to be the line of intersection of the orbital planes of the third body and the binary. Then the ascending node 
is the part of this line at which the third body crosses the plane of the binary in the direction of the angular momentum vector 
of the binary. Let £1 be the longitude of the ascending node, measured in the plane of motion of the binary (in the sense of 
the motion of its components) from the direction of pericentre of the binary. Let i be the inclination of the two orbital planes, 
and let ui be the longitude of pericentre of the third body, measured in its plane of motion from the ascending node, in the 
sense of its motion round the binary. Then the result is that the change in eccentricity is (cf. Appendix Al, eq. [A4]) 



15 m 3 / h V' 2 e\/l — e 2 
TVM 12 M 123 UJ (l + e') 3/2 



x < sin 2 i sin 20, 



.... ^ i l* _ -n3/2 

(-1/e') + \/e' 2 - 1 + - [(l + cos 2 i)cos2wsin2fi + 2cosisin2a;cos2fi] 



(7) 
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This result is admittedly quite inelegant, but the result for parabolic motion of the third body is neater; in that case 

a e = _jg f 2m y Y /2 e yr^2- sm2 nsin 2 ». ( 8 ) 

16 \M 12 3M 12 rl J v K 1 

The stated dependence on the ratio of the distances a/r p is easily understood by examination of equations (4) and (6). 
The quadrupole term in equation (4) is proportional to a/r p , and so the terms in equation (6) are proportional to a? /(Tr p ), 
where T is the period of the binary. By Kepler's Third Law the duration of the encounter is of order T(r p /a) 3 ^ 2 , and so the 
integration of equation (6) with respect to time gives a result proportional to (a/r p ) 3 ^ 2 , as we have found. 

The dependence on in equation (8) is also not hard to understand, at least for equal masses. Consider the case ui = 0, 
i — 90°, for example, and imagine the orbit-averaged binary as a distribution of mass extended along its major axis. If fi = 
or 90° then the orbit of the third body is situated on a plane of symmetry of the binary, and can exert no net torque. This 
already suggests a dependence on sin2fi. Furthermore, if = 45°, the torque increases the angular momentum of the binary, 
decreasing its eccentricity, which is consistent with the sign of the right side of equation (8). On the other hand neither the 
dependence on i nor the independence on ui can be understood from such considerations, as equation (7) shows that the 
situation for hyperbolic encounters is different. 



2.3 Initially Circular Binaries 

2.3.1 Secular Octupole Perturbation: The Power-Law Regime 

Unfortunately this lowest-order result vanishes if the initial eccentricity is zero. This is evident from equation (7) and can 
also be understood from a rather general point of view, as discussed in Section 2.3.3 and Appendix B. For a non-trivial result 
in this case it is necessary to include the octupole contribution to the interaction, i.e., the term with n = 3 in equation (4). 
Details of the calculation are given in Appendix A2, and the result for the eccentricity induced in an initially circular binary 
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may be stated as follows: 

Se _ 15 m z \m 1 - m 2 \ / M 12 \ 
VM123/ 



1/2 
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4 15 



/i(e')(!- jsin 2 i) + ^(e' 2 



1)°^ 2 (1 — 5 sin 2 oj sin 2 i) 



+ 
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where 



/i(e') = e' 4 arccos(-l/e') + ^^^(-2 + 9e ' 2 + 8e ' 4 )- 

15 

Again the result for parabolic motion of the third body is less inelegant. In fact for this case we have 



Se = 



15-7T m-Ami 



32 



m 2 \ ( Mi 2 \ 
I2M123J 



1/2 



5/2 



2 • • 2 15 . 2 -\2 2 /-1 5 . 2 -\ 
cos % sin u(l — sin 1) + cos w(l — - sin «) 



1/2 



(9) 



(10) 



(11) 



As before, it is easy to understand the dependence of these results on a/r p . The absence of O in these results should also 
come as no surprise, since the binary is assumed to be circular, and we have averaged over the position of its components. 
Another point of interest is the fact that the results are proportional to \m\ — m 2 |, and the reason for this is easily seen by 
examination of the n = 3 (octupole) term in equation (4). 

In general, for a circular binary with unequal masses, the result (9) is the dominant contribution to the induced eccentricity 
at sufficiently large r p . We refer to this as the "power-law regime", since the change in eccentricity is still proportional to a 
power of r p /a, although it is a higher power than in the case of non-zero initial eccentricity. In contrast, at smaller r p , or in 
the special case of equal masses, we show in the next section that the induced eccentricity drops exponentially with r p . 



2.3.2 Non-Secular Quadrupole Perturbation: The Exponential Regime 



In Section 2.2 we showed that, if we average over the motion of the binary and take only the quadrupole interaction, then 
the eccentricity induced in an initially circular binary vanishes. At small r p , however, the assumption on which the method of 
averaging is based breaks down, and here we re-examine the problem without averaging. The result turns out to be important 
in applications. The derivation is deferred to Appendix A3, where it is shown that 
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Here O is defined somewhat differently than in Section 2.2: it is measured in the plane of motion of the binary, from its 
position at the time of closest approach of the third body. Therefore this definition, and this formula, includes a dependence 
on the phase of the binary. Again this result simplifies a little for parabolic motion of the third body, to 
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(13) 



The approximations used in the derivation of these formulae again make use of the assumption that r p ^> a. Indeed for 
sufficiently large r p these results are negligible in comparison with equations (7) and (8) if the binary is initially non-circular, 
or with equations (9) and (11) if e = initially. Nevertheless it turns out to be an important pair of results, which are 
dominant in a small but significant range of near-encounter distances where the exponential term is not so small as to be 
masked by the previous results. 

There is one situation in which equations (12) and (13) appear to dominate at all r p 3> a: the case of a circular binary 
with equal masses. This is not a rigorous result, but we do show in Appendix B that the induced eccentricity vanishes with 
increasing r p more quickly than any power of a/r p . We already remarked that the octupole results, equations (9) and (11), 
are proportional to the mass difference |mi — m 2 \. 

Looking at equations (12) and (13) one is tempted to conclude that the eccentricity is a kind of adiabatic invariant. It 
is usual, however, to apply this term to a dynamical quantity whose variation is small when the system is subject to a large 
but slow perturbation. Here the perturbation is not only slow but small. Furthermore, the exponential form of equations (12) 
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and (13) only applies to a tiny set of orbits in which the initial eccentricity vanishes, whereas adiabatic invariants enjoy 
their special properties more globally. Henceforth we refer to this result simply as "the exponential formula" , and when it is 
dominant we say that we are in the "exponential regime" . 

It is not as easy as in the previous cases to understand the dependence of these results on a/r p , except for the exponential 
term. The exponent is simply of order the ratio of the two main timescales in the problem: the period of the binary and the 
duration of the encounter. Such a term arises in other well studied problems of this kind (cf. Goldstein 1980, eq. [11-138]). 

Finally, although the results in equations (12) and (13) vanish in the retrograde case (i — 180°), this simply means that 
the leading term in this case is of a higher order in a/r p , and we have neglected all such terms in the derivation of these 
results. 

2.3.3 Remarks on First-Order versus Second- Order Perturbations 

Details of the derivation of the foregoing results are given in Appendix A, but it is worth exposing here one or two analytical 
issues which help to clarify the nature of the results. 

There are several methods by which the computation of the change in eccentricity might be attempted. As already stated 
in Section 2.1, we have used the eccentric axis (Pollard 1976). Other possible methods are the use of Lagrange's planetary 
equations, though in their usual form they are unsuited for computations of low-eccentricity orbits, and it is better to express 
them in terms of Poincare variables (Plummer 1918). From a physical point of view the most obvious method is to compute 
the change in angular momentum. This involves the changes in both e and a, but the latter is already known (Heggie 1975, 
eq. [5.40]), and so the change in e alone may be determined. 

Whichever approach is adopted, the next step is a method of successive approximations. At lowest (first) order, unper- 
turbed Keplerian motions are substituted into the perturbation terms. At next order, the first-order results are substituted, 
and so on. In addition, it is possible to express the perturbation term itself as a multipole expansion. 

As we have seen, if we work to lowest order, retain only the quadrupole interaction, and average over the fast motion of 
the binary, then the eccentricity induced in a circular binary vanishes. It is tempting to offer the following explanation. After 
averaging, the binary components are effectively replaced by two rings of matter, which coincide if mi = mi. The density of 
each ring is proportional to the time spent by a component at each point. For a circular binary this ring is uniform and circular, 
and so the force field experienced by the third body is axisymmctric. Therefore the axial component of angular momentum 
of the third body (i.e., the component perpendicular to the plane of motion of the binary components) is conserved, and 
it follows, from overall angular momentum conservation, that the axial component of angular momentum of the binary is 
also conserved. It follows that the change in the magnitude of the angular momentum of the binary also vanishes (to first 
order). Since the change in semi-major axis over the entire encounter is exponentially small (Heggie 1975), it follows that the 
change in eccentricity also vanishes to first order. Note that this argument applies even without an expansion in powers of 
a/r p ; it depends only on the use of averaging and the use of first-order perturbation theory (whereby it is assumed that the 
interaction is calculated as if the three bodies proceed on their unperturbed orbits.) Therefore it also follows that the octupole 
contribution would vanish to first order, and indeed all terms of the multipole expansion. 

At second order it is not so clear what to expect. Averaging at first order ignores the fact that the third body causes 
the eccentricity of the binary to oscillate, and when this is taken into account (at second order) it is unclear whether the 
averaged mass distribution of the components is still axisymmetric. In their interpretation of their numerical results, Hut 
&: Paczyhski (1984) considered that the change in eccentricity must indeed be a second-order effect. We have computed the 
second-order change in the eccentricity (via the change in the angular momentum), and find a non-zero result in general in 
the case mi =fc mi. This result is generated by the octupole term. 

It turns out, however, that the distinction between first- and second-order effects is not an absolute one, but depends on 
the variables used in the computation. In the present paper we adopt the Lenz vector in preference to the angular momentum, 
as a means of computing the eccentricity. With this variable it turns out that a non-trivial result for the eccentricity induced 
in a circular binary is obtained at first order (provided again that mi 7^ mi.) The result agrees with that obtained from 
a second-order computation of the change in angular momentum, and is also caused by the octupole interaction, but it is 
much more easily calculated using the Lenz vector, as only first-order perturbation theory is needed. It is not hard to see why 
second-order perturbation theory is required if the angular momentum is used. The angular momentum is proportional to 
\/l — c 2 , and so if the initial eccentricity is zero, the change in angular momentum is second order in the change of eccentricity. 

Finally in this section we turn to the case of a binary with equal masses. Here the odd-order terms in the multipole 
expansion, equation (4), vanish, and it is easy to see that all first-order calculations, whether with the Lenz vector or the 
angular momentum, will yield a null result. Consider, for example, the method adopted in Appendix A. The perturbation is 
an odd function of r (the position vector of one component relative to the other), and it follows from equation (6) for the rate 
of change of the eccentric axis that this is also odd. A null result is obtained, therefore, if we average over the motion of the 
components. This result, which is so simple to prove at first order, also applies at all orders of perturbation theory, but the 
proof is somewhat technical, and is relegated to Appendix B. The implication of the result, however, is that the change in 
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eccentricity must vanish more quickly (with increasing r p ) than any power of the expansion parameter a/r p , and it is likely 
that the result vanishes exponentially as an inverse power of this parameter (cf. Section 2.3.2). 

The problem of rigorously justifying the computation of exponentially small perturbations is one which has attracted 
considerable attention from mathematicians in recent years (e.g., Byatt-Smith & Davie 1990). The difficulty is, of course, 
that standard perturbation methods proceed by expanding in powers of a small quantity, whereas these exponentially small 
results have no power series expansion. Unfortunately it still appears that each case must be treated separately if it is to be 
handled rigorously, and there is no general calculus for this class of problems. The present problem is certainly considerably 
more elaborate than those which have been treated properly. Our interest in this result is more pragmatic, however. We shall 
see that the result of Section 2.3.2, which exhibits exponential dependence, is important when it dominates the result from 
the octupole interaction (Section 2.3.1), and in this situation it does not require elaborate justification. 

2.4 Illustrations and Comparison with Numerical Results 

We have performed a large number of numerical integrations of encounters between a binary and a third star. In this section 
we present a sample of results, to establish the accuracy of our analytical expressions and to illustrate the results for some 
typical cases representative of the different possible regimes. 

We shall focus primarily on parabolic encounters, which are particularly relevant astrophysically since the conditions in 
globular clusters are very close to the parabolic limit except for very wide binaries (a> 1 AU). In general, in any star cluster, 
hyperbolic encounters are of importance only for binaries that are sufficiently wide to be classified as "soft" and therefore 
will tend to be easily disrupted by interactions (Hills 1975; Heggie 1975; Hut 1983). Deviations from the parabolic results 
remain small as long as the relative velocity at infinity for the encounter is less than the orbital velocity of the binary, and 
the encounter is not too distant (cf. eq. [1] and Section 3.2 below). 

2.4-1 Numerical Integrations 

The numerical integrations are fairly straightforward. We have used the Bulirsch-Stoer algorithm described in Press et al. 
(1992) to integrate the differential equations of motion of the three-body problem directly. Total energy and total angular 
momentum conservation is maintained typically to within 10~ 12 — 1CP 11 , for computations in double precision. This allows 
changes in eccentricity as small as 8e ~ 10 -8 to be measured accurately. The input parameters for a given integration are 
the masses mi, m,2 and mz, the semi- major axis a, eccentricity e, and mean anomaly <j)i (or orbital phase, measured from 
pericentre) of the binary at t = 0, the periastron distance r p and velocity at infinity V characterizing the outer orbit (cf. 
Section 2.1), the initial separation n (distance between mz and the centre of mass of the binary at t — 0), and the angles 
£7, i, and lo specifying the relative orientation of the two orbits (defined in Section 2.2). We generally adopt units such that 
G = mz — a — 1. Each integration follows the system through periastron passage and is continued until the final separation 
has increased to a value equal to the initial separation, r/ = r;. Typically we use ri/r p ~ 10 2 . For a given set of input 
parameter values, we normally repeat the integrations for many different values of the initial phase 0;, covering systematically 
the interval from to 2n. This allows us to compute a phase-averaged result and to study the extent of the phase dependence. 

2.4-2 Dependence on the Periastron Distance 

First we consider three simple situations where one of the three expressions (8), (11), and (13) dominates over the entire range 
of periastron separations where 8e remains small enough that a perturbative calculation applies. 

In Fig. 1, we show the change in eccentricity of a binary with non-zero initial eccentricity, following an encounter with a 
star on a general inclined orbit. The three masses are assumed identical. The result is given by equation (8) to high accuracy. 
Here and in all subsequent figures the solid lines show the analytical results and the dots show the numerical results, with 
the "error bars" indicating the full extent of the dependence on orbital phase. In this case the phase dependence remains 
completely negligible down to very small values of r p /a. Notice that, as long as the three masses are comparable, the magnitude 
of the fractional change in eccentricity \5e\/e predicted by equation (8) is always very small. Using equation (8) we see that 
a change |<5e|/e <; 1 would require mz <; (r p /a) 3 Mi2. 

Now we turn to initially circular binaries. To isolate the octupole effect (equation [11]), we consider a binary with mi 7^ 7712 
and we focus on distant encounters, with r p /a <; 5. The variation of the induced eccentricity 5e with r p is shown in Fig. 2 for 
a coplanar, prograde encounter, and for two different values of m,2. Notice the steeper power-law behavior (r p °^ 2 ) compared 
to the case of an initially non-circular binary (r p 3 ^ 2 ). Notice also the large drop in Se when mi becomes nearly equal to mi. 
To best illustrate the exponential regime (equation [13]), we show in Fig. 3 the variation of Se in the equal-mass case. The 
induced eccentricity now drops much more rapidly, to values <C 1CF 7 for r p /a J> 10. Notice that equation (13) continues to 
predict Se with remarkable accuracy down to very small values of r p /a. 

In a general case, all three types of perturbations dominate in different ranges of r p . Consider for example a typical 
binary millisecond pulsar containing a neutron star of mass mi = 1.4 Mq with a low-mass companion mi = 0.2 Mq in a 
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Figure 1. Change in eccentricity 8e for an initially non-circular binary (initial eccentricity e = 0.1 or 0.01) perturbed by an encounter 
with a third star on an inclined, parabolic orbit with periastron distance r p (given in units of the binary semi-major axis a). The dots 
show the results of numerical integrations; the lines show the analytical power-law result, equation (8). The "error bars" indicate the full 
extent of the phase dependence. For r p /a <; 5 this phase dependence is no longer visible. 

nearly circular orbit of initial eccentricity e = 10~ 6 ; the binary is in the dense core of a globular cluster where the average 
stellar mass is = 1 Mq (cf. Rasio & Heggie 1995). Fig. 4 shows the final eccentricity e/ = e + Se following a coplanar, 
prograde encounter, as a function of r p /a. The exponential regime corresponds to r p /a <! 5. For r p /a <; 200, equation (8) 
dominates and the fractional change Se/e is very small. In a wide range of intermediate values, 5 S5 r p /a 200, the octupole 
power-law result is dominant. As discussed in Rasio & Heggie (1995), an important consequence of these results is that the 
cross section for inducing eccentricity in binary millisecond pulsars in clusters is much larger than was estimated in previous 
studies (see also section 3 below). 

2.4-3 Angular and Phase Dependences 

Although we have checked very carefully most of the phase and angular dependences appearing in our analytical expressions, 
we shall only present a small number of illustrative cases here. We focus on initially circular binaries, since the angular 
dependent factors appearing in equation (8), for the case of non-zero initial eccentricity, can be understood from elementary 
arguments, at least in part. We use the same typical values of the masses for a binary millisecond pulsar in a globular cluster 
as in Fig. 4, and we consider only parabolic encounters. 

First, we look at the dependence on inclination for the induced eccentricity in the power-law regime. Fig. 5 shows 
the results for r p /a = 30 and for two particular values of u>. For u> = 0, the first term in the angular dependent factor of 
equation (11) vanishes and the dependence on inclination is oc |1 — (5/4) sin i\, giving zeros of Se at i ~ 0.3527T and i ~ 0.6487T. 
The maximum Se is obtained for i — and i — tt, with a secondary maximum at i = 7r/2. For ui = 7r/2, the second term in 
the angular dependent factor vanishes and Se oc | cosi [1 — (15/4) sin 2 i\ |, giving zeros at i = tt/2, i ~ 0.1737T, and i ~ 0.8277T. 
Again the maximum Se is obtained for i — and i = ir, with two secondary maxima at i ~ 0.3357T and i ~ 0.6657T. The 
agreement between the numerical and analytical results is everywhere excellent, with deviations remaining always ^ 3%. 

In the exponential regime, a much stronger phase-dependence is observed, and larger deviations are observed between 
the numerical and analytical results. This is not surprising, since the exponential regime corresponds to rather large values of 
the perturbation expansion parameter a/r p . In Fig. 6, we illustrate the strong phase-dependence as well as the dependence on 
inclination for r p /a = 3.5. The phase-dependence is largest for i = (coplanar prograde encounters), which also corresponds 
to the maximum value of Se according to equation (13). For i = tt, equation (13) predicts Se = even though the numerical 
result is clearly non-zero. The small discrepancies between numerical and analytical results in this regime come from neglecting 
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Figure 2. Change in eccentricity Se for a binary with initial e = perturbed by an encounter with a third star on a coplanar, prograde, 
parabolic orbit with periastron distance r p (given in units of the binary semi- major axis a). Conventions are as in Fig. 1. The lines 
show the analytical result, equation (11), for two different values of ni2. The dashed line shows, for comparison, the variation of Se for a 
non-circular binary with initial e = 0.01 (same case as in Fig. 1). 




Figure 3. Change in eccentricity Se for a binary with initial e = and containing two identical masses mi = m,2- Conventions arc as in 
Figs. 1 and 2. The line shows the analytical result, equation (13). 
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Figure 4. Final eccentricity ef = e + Se for a binary with initial e = 10~ 6 following a coplanar, prograde encounter. Conventions are as 
in Fig. 1. The masses have values representative of a typical binary millisecond pulsar in a globular cluster. The solid lines show the three 
separate analytical results, equations (8), (11), and (13). The dashed line shows the continuation of the exponential result, equation (13), 
for comparison. Notice that the octupole power-law result, equation (11), dominates over a wide range of intermediate values of r p . 

the more rapidly oscillating term in the derivation of equation (13) (cf. Appendix A3), and from the octupole perturbations. 

Other circumstances where a strong phase- dependence is observed are when two different perturbation effects contribute 
similar-amplitude changes in the eccentricity, as in Fig. 4 near r p /a ~ 5, or when \5e\ ~ e. As an example of the latter 
consider the case of a binary with mi = m 2 and i = (so that both power-law expressions [8] and [11] vanish) and with a 
very small initial eccentricity e = 1CF 4 , as illustrated in Fig. 7. The principal contribution to 5e is given by equation (13), 
though this was calculated for the case of an initially circular orbit. Since in the present case e is small but non-zero, 8e must 
be interpreted as the magnitude of the vector 5b, and the final eccentricity is given by e/ = |e + Se\. Now the direction of 8e 
depends on the phase of the binary, as can be seen by careful inspection of eq.(A21), and so if |<5e| = e we see that e/ = 
for certain phases. In fact eq.(13) shows that the condition \8e\ — e is satisfied when r v /a ~ 6.2. At this value the numerical 
result in Fig. 7 would have an "error bar" extending to arbitrarily small values of e/ in this logarithmic plot. An alternative 
argument may be constructed by considering the time-reversal of such an encounter, i.e. one in which an initially circular 
binary is perturbed into an orbit with final eccentricity e/ = 1CF 4 . Eq.(13) shows that this can only occur at the same critical 
value of r p /a. 

2.4-4 Dependence on the Masses 

In general, when the three masses are all of comparable magnitude, and mi =fc mi, the mass-dependence of the induced 
eccentricity Se is rather weak. Near mi = mi, however, the power-law expression, equation (11), can become arbitrarily small. 
The peculiar dependence of 8e on |mi — m 2 \ in the power-law regime is illustrated in Fig. 8 for a typical case. The agreement 
between analytical and numerical results is again excellent. 

When m 3 M 12 , as in the case of an encounter between a binary and a massive black hole, it is possible to find that the 
exponential result dominates even at rather large values of r p /a. This is because the duration of the encounter decreases as 
m,3 increases at fixed r p . In Fig. 9 we show the dependence of 8e on m 3 for fixed r p /a = 30. The exponential result becomes 
dominant when m$ J> 300. For ms J> 3000, the binary is disrupted by the encounter. 

2-4-5 Hyperbolic Encounters 

For fixed r p , as the relative velocity at infinity V increases from zero, the duration of the encounter becomes shorter and 
shorter and we always expect the exponential result, equation (12), to dominate eventually. An example is shown in Fig. 10, 
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Figure 5. Dependence of induced eccentricity Se on inclination for an initially circular binary with mi = 1.4, rri2 = 0.2, 7713 = 1, for a 
parabolic encounter with r p /a = 30 (power-law regime). The round dots are for ui = 0, the square dots for u> = w/2. The lines show the 
analytical results, equation (11). 




Figure 6. Same as Fig. 5 but for r p /a = 3.5 (exponential regime). The two lines correspond to cos(4u; + 2Q) = ±1 in equation (13), 
i.e., they show the extent of the theoretically predicted phase dependence. As before, the "error bars" show the full extent of the phase 
dependence in the numerical results while the dots show the phase-averaged numerical results. 
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2 3 4 56789 10 

Figure 7. Same as Fig. 3 but for a binary with non-zero initial e = 10~ 4 . Notice the strong phase dependence near r p /a ~ 6 — 7, where 
the magnitude of 5e (given by eq. [13]) is comparable with the initial value of e. 




Figure 8. Dependence of Se on the masses near mi = m,2 for the power-law regime. Equation (11) predicts <5e = for mi = ni2, in 
agreement with the numerical results. 

where we fix r p /a — 30. Although the secular octupole effect is completely dominant in the parabolic limit, non-secular effects 
become dominant for V/(Gm 3 /a) 1/2 > 4, or V > 120 kms" 1 (m 3 /l M ) 1/2 (a/l AU)~ 1/2 . Between V = and this value, <5e 
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Figure 9. Same as Fig. 4 but showing the dependence of Se on mz- For sufficiently large 7713 (here for m.3 <; 300), non-secular effects 
can always become dominant, even for encounters with large r p /a. 

decreases by almost an order of magnitude, but then increases again as V increases in the exponential regime. In contrast, 
Fig. 11 shows a case where the exponential result dominates everywhere, which typically happens at sufficiently small r p . 
Here Se increases at first, reaches a maximum at V j (Gmzj a) 1 ^ 2 ~ 2 and then decreases with increasing V. Note that the 
impulsive limit, where Se oc V^ 1 (cf. Section 3.2 and Appendix A4), would only be reached for much higher values of V, and 
therefore does not seem astrophysically relevant. 



3 CROSS SECTIONS FOR ECCENTRICITY CHANGE 
3.1 Analytical Results 

While the results of Section 2 are fully detailed, giving the dependence on all the circumstances of the encounter, for the 
purpose of applications it is necessary to average over several parameters which are normally only known statistically. In the 
present section we subject the results to some further processing in order to extract the cross section for encounters leading 
to changes in eccentricity above a given value. The cross section is essentially the square of the impact parameter p, but we 
shall restrict attention to the case of near-parabolic encounters. In this case the relation between p and the distance of closest 
approach reduces to 

v 2 V 2 

^^AW (14) 
3.1.1 Binaries with Non-Zero Initial Eccentricity 

Consider first binaries with non-zero initial eccentricity. From equation (14) it follows that the dependence of equation (8) on 
the circumstances of the encounter may be summarised as 

Se — — Ap~ 3 sin2f2sin 2 i, (15) 



where 

Vl-e 2 . (16) 



15tt m 3 M 123 (Ga) 3/2 



4 s/M^ V 3 

Now we compute the cross section for events in which Se > 8e , where Seo is some value of interest. We shall suppose that 
8e > 0, and then it will be obvious how to extend the result to negative values, as expression (8) obviously takes positive 
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Figure 10. Variation of the induced eccentricity <5e as a function of V, the relative velocity at infinity for an encounter with fixed 
r p /a = 30. The two solid lines show the power-law and exponential analytical results (equations [9] and [12]); the dashed line shows their 
sum. The exponential result dominates for VKGmz/a) 1 / 2 <; 4. 




Figure 11. Same as Fig. 10 but for r p /a = 4. The exponential result (upper solid line) dominates everywhere in this case. Notice how 
the variation of <5e is opposite to that shown in Fig. 10, with 5e increasing at first, reaching a maximum for V / {Gm^,/ a) 1 / 2 ~ 2, then 
decreasing at higher velocities. 
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Table 1. Correction to the Asymptotic Cross Section in the Exponential Regime 
Mio3<5eo 

ll 0.333 0.25 0.2 0.1 0.033 0.01 0.001 0.0001 0.00001 

X 0.112 0.228 0.307 0.497 0.705 0.850 0.938 0.968 0.982 

and negative values with equal probability. Using the appropriate statistical distributions of i and Q, we write 

a(5e > Seo) = / 2%pdp—dQ- sin idi, (17) 
J 2n 2 

where the domain of integration is given by p 3 < — A sin 2Q sin 2 i, < i < tt and, for fi, the union of the intervals [ir/2, 7r] and 
[37r/2, 2ir] (to ensure that Se has the same sign as 8eo). The integrations are straightforward, and lead to the result 



a(Se > Seo) = 



-<!Mf 



2 



(18) 



V Seo J 147r 

After the definition of A is reinstated the result becomes 
/x . . V3 /15tt\ 2 / 3 r /2\ /5\1 2 / ml \ 1/3 GM 123 a 2/3 2,1/3,-2/3 

It is clear that the cross section for negative changes, i.e. <5e < — Seo, where Seo > 0, is given by the same expression. 

3.1.2 Initially Circular Binaries 

Now we turn to the cross section corresponding to the parabolic result for a circular binary in the power-law regime, equa- 
tion (11). Everything is as straightforward as in the previous case, with some minor alterations. First, only positive values 
of <5eo are relevant. Second, the appropriate distribution for uj is uniform on the range [0, 2ir]. Third, the integration cannot 
apparently be performed analytically, and it is a tricky job numerically because the integrand is not smooth at certain points 
in the relevant domain. For this reason we have opted for a Monte Carlo evaluation, and so we find that 

a(Se > Se () ) * 4.62 ^ j j 5e . (20) 

The coefficient is correct to the stated number of figures. 

The exponential regime for circular binaries is somewhat trickier, because the dependence on r p in equation (13) is not 
a simple power law. Using equation (14) we first recast equation (13) as 

Se = A(Bp 3 ) 1/2 exp(-Bp 3 )/2(i, u, Q), (21) 



where 



A = 3V6WMi2 3 , ^4£) 1/2 (GmLV 2 ' (22) 

and /2 is the geometric factor at the end of equation (13). Proceeding as in the previous cases, we can write the desired cross 
section as 

a(Se>Se )= / 2npdp^- smidi-^-dO, (23) 
J 2 2tt 

where we have set 9 — 4u) + 2Q: it is clear that the appropriate distribution is uniform. The integration with respect to p gives 
a factor J pdp = p 2 /2, where p is the root of equation (21) with Se replaced by Seo- It follows that 

a(Se > Se ) = J x 2/3 smidid6, (24) 

where x is the larger root of x — (1/2) lna: = \n{Af2/8eo). An approximate solution is x = \n(A/Seo), which is valid when 
ln(^4/5eo) ^> 1, and in that limit we have a(Se > Seo) = 7r[B _1 ln(A/Seo)] 2 ^ 3 , i.e. 

a{S e>Seo)=^ — ) -^-^^-^-^ x (_j. (25) 

Here \ is a correction factor which will be used to allow for a range of values of Seo /A; in the limit of Seo «lwe use x(0) = 
1. For other values of Seo /A we have evaluated the integral in equation (24) by a Monte Carlo technique, with results which 
are given in Table 1 as the correction \ to the foregoing asymptotic result. The data are accurate to 1 unit in the last digit 
printed. 

The maximum possible value of Seo is A/\/2e, ~ 5.6m 3 /Mi23- 
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3.2 Comparison with Previous Work 

It is interesting to compare these cross sections with those published by other authors. 
Hut & Paczyriski (1984) gave a result which may be expressed, using our notation, as 

IS ^ X \ I a v _i_ n n\ GM l2 3a ^mzM 12 _i/ 3 

o H p(Se > Se ) = (4.7 ±0.9) — — M Se Q . (26) 

This was obtained by numerical scattering experiments for hard, initially circular binaries, and so we shall compare it with 
our equation (20). The result is 

^ = (1.02 ± 0.19) , mrM }i, l)H Sel /15 . (27) 

According to Hut & Paczyriski (1984), their result is valid to within the stated accuracy of 20% for certain ranges of a and V, 
and for mass ratios in the ranges 0.03 < mi/mi < 0.3, 0.2 < mz/mi < 1. In this range of masses, the mass-dependent factor 
in equation (27) varies steadily from about 0.59 (at mi = 0.03mi, m 3 = mi) to about 0.95 (at mi = 0.3mi, mz = 0.2mi). 
Actually for the most discrepant result here the scattering cross section of Hut & Paczyriski applies to binaries which include 
the borderline between hard and soft pairs, because m 2 is so low. Therefore it is probable that their cross section is a 
satisfactory compromise in this domain. Certainly, our own result is not intended to be applicable to such encounters. 

Now we turn to a detailed comparison with a representative result from the paper by Rappaport et al. (1989, hereafter 
RPV). They give results for a single-parameter family of circular binaries in each of two environments, characterised by 
mixtures of stars with given mass, number density and velocity dispersion. We have chosen to compare with their result for 
the hardest binary (one of orbital period 3 days) in their model for the stellar population of the star cluster co Cen. Our 
first step is to average our result, equation (20), over a Maxwellian distribution of velocities, though with an extra factor 
of V in the averaging to account for the enhanced rate of interaction between stars of high relative speed (as was done by 
RPV). The result is that the factor V 2 in equation (20) is replaced by 2({v 2 ) + (w|))/3, where the averages are the mean 
square three-dimensional speeds of the binary and the third body, respectively. Next the cross section is averaged over the 
stellar species in the mixture specified by RPV. Then a similar process is carried out for the cross section, equation (25), 
corresponding to the exponential regime. This is slightly more awkward, as the averaging over the stellar species must be 
carried out afresh for each value of Seo- For comparison with the results of RPV we simply adopt the larger of our two cross 
sections at each value of Seo- 

Fig. 12 shows the ratio of RPV's cross section to ours as a function of Seo- The change in slope near Seo m 0.003 is 
associated with the point at which the cross section of equation (25) begins to dominate that of equation (20). The change at 
Seo = 0.01 has a similar origin in the work of RPV: the functional form of their cross section changes at this point. From there 
up to the point where Seo — 0.5 our cross section is systematically higher, by about 60%, than that of RPV. The reason for 
this systematic disagreement is not certain, but it is clear that our approach of using the maximum cross section at each value 
of Seo is rather crude. The sharp drop towards as Seo — > 1 reflects the fact that our cross section does not vanish in this 
limit, whereas that of RPV does so. Turning now to very small values of Seo, we see that the cross section of RPV decreases 
monotonically relative to ours. Indeed, their expressions are similar to our equation (25), valid in the exponential regime but 
not in the power-law regime. This may be related to their use of relatively short and low-precision numerical integrations, 
which are probably inadequate for computing the very small changes in eccentricity corresponding to the power-law regime. 
It must also be pointed out, however, that our cross sections were calculated for parabolic encounters, which is inappropriate 
when r v /a J; GMnz/iaV 2 ). Substitution of typical values into equation (11) therefore shows that the validity of our result is 
restricted to the range Seo iJ 10~ 5 in this case. 

Next we turn to the work of Hills (1991), who computed cross sections for equal-mass binaries in encounters with more 
massive intruders (mz mi = mi). The initial eccentricity was e = 10~ 4 , and since Hills' results correspond to much 
larger changes, in the range Se > 0.1, we shall compare equation (25) with his results extrapolated to V = 0. He gives a 
scaled cross section a(8e > 0.3) for mass ratios in the range 10 < mz/m\ < 10 4 . Our results show the same trend with 
increasing mass ratio, but are consistently smaller by about 10-20%. The dependence of the cross section on Seo can also be 
checked, as Hills gives the ratios of a(8e > Seo) for the values Seo = 0.1, 0.3 and 0.7. For the smaller values of Seo our result 
agrees with Hills' to within the same accuracy of 10-20%. For the larger values the agreement is poorer: Hills finds values 
of a(Se > 0.3)/a(5e > 0.7) in the range 1.5-3, whereas our results yield a ratio closer to 1.4. Since our results are based on 
assumptions which are valid for small eccentricity changes, it is perhaps surprising that the agreement is so good. 

Finally we compare our results with those of Phinney (1992), who gives formulae for the induced eccentricity in an 
initially circular binary, and an approximate cross section. Like us, Phinney assumes that r p 3> a. For simplicity, he only 
considers the coplanar case and he assumes a straight-line (i.e., extremely hyperbolic) trajectory. In this limiting case, the 
exponential result always dominates and the appropriate comparison is to our equation (12) in the limit where e' — > oo and 
for i = 0. Identifying Phinney's quantity A with (Mii/M 1 i 3 ) 1 ^ 2 (r p /a) 3 ^ 2 (e')~ 1 ^ 2 in our notations, we find perfect agreement 
with Phinney's equation (5.2). Note that Phinney's equations (5.2) and (5.3) have been mislabelled: they correspond to the 
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Figure 12. Comparison between the scattering cross sections cr(8e > 5eo). The result of Rappaport et aJ. (1989), denoted by the 
subscript RPV, is compared with the cross section derived from the results of this paper, for a range of values of 5eo- 

prograde and retrograde cases, respectively (Phinney, personal communication). Phinney's result for the coplanar retrograde 
case is of higher order and was neglected in our derivation of equation (12). 

Phinney (1992) also gives a result for fast encounters (the "impulsive limit"), which we now consider briefly, although we 
do not expect this case to be of importance in most astrophysical applications (one possible exception being the perturbation 
of wide circular binaries by halo stars in the Galaxy). We give a derivation in Appendix A4, and our result for coplanar 
encounters with an initially circular binary is 



where 9 is the angle between the vectors r and R at the time of pericentre passage. This result vanishes when 9 — 90° , which 
can be understood because in that case the third body generates equal impulses in both components, and therefore does not 
alter their relative motion. Phinney's equation (5.1) as published differs in the trigonometrical factor. This is because of a 
typographical error: the first term inside the square root, which is printed as 3sin 2 #o should have been 3sin 2 2#o (Phinney, 
personal communication). With this typo corrected, it is easy to show that Phinney's equation (5.1) is in fact identical to our 
equation (28). 

4 CONCLUSIONS AND DISCUSSION 

We now summarise our main results. They describe the change Se in the eccentricity of a binary of semi-major axis a which 
experiences an encounter with a third body, in the following regime: if r p is the distance of closest approach of the perturber, 
we suppose that r p considerably exceeds a and that the time scale of the encounter considerably exceeds the period of the 
binary. With these assumptions our results are of most relevance to hard binaries, and the orbit of the third body is nearly 
parabolic. 

1. When the induced change in eccentricity, Se, is small compared with the initial eccentricity of the binary, then Se falls off 
as (a/r p ) 3//2 . The relevant results are equation (7) for a hyperbolic encounter and equation (8) in the parabolic limit. Other 
factors affecting the result are the masses of the participants, the geometry of the encounter, the initial eccentricity of the 
binary and the orbit of the perturber. The result comes from the quadrupole interaction between the binary and the third 
body. Its main relevance is to binaries with a large initial eccentricity (Fig. 1). 

2. When Se is large compared with the initial eccentricity, and r p much exceeds a, the final eccentricity falls off in general 
as (a/r p ) 5 / 2 (equations [9] and [11]), and also depends on various other factors. The fact that this result originates in the 




(28) 
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octupole interaction accounts for the extra power. This result is of relevance mainly to binaries whose initial eccentricity is 
small (e.g., initially circular binaries), except for very close encounters (cf. Figs. 2 and 4). 

3. When r p exceeds a, but not by a large factor, the induced eccentricity falls off exponentially with a power of r p /a, though 
the coefficient also depends on this ratio, and other factors (equations [12] and [13]). The exponential factor occurs because, 
at these values of r p /a, the time scale of the encounter begins to be comparable with the period of the binary. It is mainly 
of relevance to binaries with small initial eccentricity, and to encounters which lead to a large induced eccentricity, but at 
periastron distances r p where the quadrupole interaction is still dominant. It is also relevant to the very special case of a 
binary with components of equal mass (Fig. 3), or to the more important case of a massive perturber (Fig.9). 

4. Cross sections corresponding to these three results are given in equations (19), (20) and (25). The first two depend on 
certain powers of the change in eccentricity, while the dependence in the third is approximately logarithmic. 

These results have been compared with numerical data in various ways. We have tested the predicted dependence of the 
induced eccentricity change Se on the encounter distance r p (Figs. 3, 4 and 7), on the masses of the stars (Figs. 8 and 9), and 
on the geometry of the encounter (Figs. 5, 6, 10 and 11). We have almost always found excellent agreement between numerical 
and analytical results. 

We have also compared our results, and especially the cross sections, with those derived by other authors. The main 
point at which previous work requires revision is the behaviour of initially circular binaries. Previously it was thought that 
the induced eccentricity 8e falls off exponentially with encounter distance, but it is now clear that this result applies only to 
a limited range of periastron distance, beyond which 8e falls off only as a power of the distance. The effect of this is that 
the cross section for encounters leading to small (but observable) induced eccentricities is very considerably larger than was 
previously sometimes thought. 

The advantage of using our formulae for the induced eccentricity change in future investigations is speed. Where the 
parabolic approximation is valid, the stated cross sections may be used directly, and it is a straightforward matter to compute 
the cross section for any reasonable mass spectrum. Where the parabolic approximation is invalid, but the encounters are still 
both slow and tidal, cross sections may be computed in a Monte Carlo fashion, the outcome of each encounter being obtained 
from the appropriate formulae of Section 2. In this way estimates of cross sections in a wide variety of collisional stellar 
systems may be estimated much more rapidly than in the past, without the need for numerous careful and time-consuming 
three-body integrations. 

Currently the most important application of our results concerns low-mass binary millisecond pulsars (LMBPs) in globular 
clusters. These are formed when an old neutron star in a binary system accretes material from a red-giant companion. This 
leads to the production of a spun-up (recycled) pulsar with a low-mass white-dwarf companion in a wide, circular orbit (typical 
eccentricities at birth being in the range e» ~ 10 -6 — 10~ 3 ). This formation process is well understood theoretically (Verbunt 
1993; Phinney & Kulkarni 1994). In globular clusters, the wide orbits of LMBPs can be perturbed significantly by passing 
stars. Thus final eccentricities e/ » e, can be observed today, and the measured values contain important information about 
the dynamical history of the binaries and their environment. In fact Rasio & Heggie (1995) have shown, using the analytical 
results derived here, that all currently known LMBPs in globular clusters must have been affected by interactions, with their 
current eccentricities being at least an order of magnitude larger than at birth. 

Distant interactions with passing stars may also be important for the orbital evolution of X-ray binaries in globular 
clusters. Hut & Paczyhski (1984) have shown that even a very small change in eccentricity can lead to a catastrophic increase 
in the mass transfer rate in a semidetached binary. In some cases, an eccentricity change as small as 8e ~ 10~ 4 — 10~ 3 
could even lead to the formation of a common envelope, ultimately destroying the binary. Other problems where our results 
could find applications include the interactions of binaries with massive black holes (Hills 1991) and the orbital evolution of 
a massive black-hole binary in a galactic nucleus (Mikkola & Valtonen 1992). These problems involve extreme mass ratios 
(m 3 3> M12 or m 3 <C M 12 ) which we did not address specifically in this paper, but our analytical expressions should apply in 
those cases as well (cf. Fig. 9) . 
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APPENDIX A: DERIVATION OF THE ANALYTICAL RESULTS 
Al Nonzero Initial Eccentricity 

Our task here is to derive equation (7) by integration of equation (6), where F is given by equation (4). At lowest order, which 
is all we shall need, we substitute unperturbed expressions for r and R on the right side of equation (4). 

We denote by a and b two orthogonal unit vectors in the original plane of motion of the binary, such that a is in the 
direction of pericentre. (In fact e = ea.) We also let c = a x b be the unit vector in the direction of the angular momentum. 
Similarly we define two vectors A and B in the original plane of motion of ms relative to the barycentre of the binary. 

Next, define coordinates x,y,X and Y such that r = xk + yh and R = A A + VB. Then expressions for these four 
coordinates in unperturbed Keplerian motion can be found in standard texts (e.g., Plummer 1918; Danby 1962). Thus 
x = a(cosE — e) and y — bsinE where b — a\/l — e 2 and E is the eccentric anomaly. It is related to time by Kepler's 
equation, n(t — to) = E — e sin E, where 

nV = GM 12 , (Al) 

and to is the time of pericentre passage. If the relative motion of the third body is hyperbolic with eccentricity e', for this 
section and the next one the expressions are most conveniently written in polar form, i.e., X — RcosO and Y = RsinO. Here 
the polar angle is related to time implicitly via angular momentum conservation, i.e. R 2 9 = [GMi23r p (e' + l)] 1 ^ 2 , where r p 
is the distance of closest approach of the third body to the barycentre of the binary. We take t = at this instant, i.e., when 
6» = 0. Also R = r p (e' + l)/{e' cos 6 + 1). 

The integration of the equation for e proceeds by first averaging over the motion of the binary (at fixed R). At this point, 
if only the quadrupole term (n = 2) is retained in equation (4), the result is 

<e> = -jg^ [15£„ a + (-9 ? 2 + 6 V 2 + 3< 2 ) b + 3^c] , (A2) 

where T is the period of the binary, and £, r\ and ( are defined implicitly by R = £a + rjh + £c. In fact only the component 
of equation (6) along a is needed, since eSe = e ■ Se = ea • 8e; and since Y is an odd function of t we find that the change in 
eccentricity is 

Se = -™E*L«^ fa-Ab-A r ^dt + a-Bb-B /°° ^dt) . (A3) 
M12 T \ J^R* J^R* J 

For hyperbolic motion of the third body the result is 



2e' 2 V Mi 2 Mi 23 rg(l + e') 3 

ja-Ab-A 3e' 2 arccos + (4e' 2 - l)Ve' 2 - 1 +a-Bb-B 3e' 2 arccos (-^7) + (2e' 2 + 1)\/V 2 - lj j 



(A4) 
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For parabolic motion (e' — > 1) the corresponding result is 
15tt / 2m 2 a 3 N 1/2 



5e = -T(i«|J «>A^(^Ab.A + & .Bb.B) > (A5) 

a result already given in Heggie (1975, eq. [5.66]), though with the wrong sign. It is convenient to write the direction cosines 
a • A, b • A, etc., in terms of the Keplerian elements outlined in the text. This leads to equations (7) and (8). 

A2 Initially Circular Binaries: The Power-Law Regime 

Both of these results vanish for an initially circular orbit. However, the octupole (n = 3) terms in equation (4) yield a non- 
trivial result in this case, which is the case we consider henceforth. As before we take t = at the instant of closest approach 
of rri3, but now we can choose a to be in the direction of r at that time. The analogue of equation (A2) is 

... 15m3 (mi - m 2 )7ra 4 / <" 2 \ , - t C\ 

<e) = — mww> — I, 1 - °jp ) {r > a - ?b) - 

Converting to coordinates in the plane of motion of the third body, we find that 

Se = f r W dt I 1 - 5 ^-) 2 ] A xc+ /°° ^fdt {[1 - 5(B . c) 2 ] A x c - 10B . c A . cB x c} ) .(A6) 

12 \J—oc J-cc / 

Note that Se lies in the plane of motion of the binary, and is independent of the choice of a and b in this plane. Incidentally, 
it is also easily seen that <5e is orthogonal to A in certain circumstances, e.g., in coplanar motion. This means that the binary 
orbit becomes elongated in a direction orthogonal to the position of 7713 at pericentre. 

For hyperbolic relative motion of the third body the two integrals in equation (A6) have the values 



X 3 _ y/e' 2 - 1(2 + lie' + 32e' ) + 45e' arccos(-l/e') 
R7 ~ 30e' 3 x/GM 123 4(e' + 1)5 



and 



f 

J OO 



XY 2 _ \Je' 2 - l(-2 + + 8e' 4 ) + 15e /4 arccos(-l/e / ) 
R7 SOe'VGMiasr^e' + l) 5 



From these results the values of 5e and Se = \5e\ are easily derived, though the expressions are relatively cumbersome. It 
helps to resolve Se along and perpendicular to the line of nodes, and to express the various products of unit vectors in terms 
of the elements i and lo already introduced in §2.2. Eventually we obtain equation (9), where /i(e') is defined in equation (10). 
Note that, in the case of parabolic motion of the third body, we have /i(l) = 7T, and so the result is rather more graceful. 

A3 Initially Circular Binaries: The Exponential Regime 

As we have seen, averaging over the fast motion of the binary leads to the conclusion that the quadrupole term gives a 
null result in the case of an initially circular binary, and that the octupole term provides the leading term in the induced 
eccentricity. At a sufficiently small distance of closest approach, however, the method of averaging becomes inaccurate, and 
so a different approach to the evaluation of the quadrupole term is needed. We continue, however, to restrict attention to an 
initially circular orbit. 

In this case the contribution from the quadrupole term is 



m 3 



M 12 R 3 



g (r-R)(r-R) (r ■ R) 2 . 2 . 
6 R* r-3-^r + 7-r 



(A7) 



As before we take t = at the time of closest approach of 7713, and choose the vector a to be the direction of r at this instant, 
so that r = a(a cos nt + bsinnt), where n is the mean motion of the binary. 

Since the right side of equation (A7) involves the third power of components of r or r, its terms involve factors of exp(±int) 
and exp(±3int) when written in exponential form. From what follows we shall see that the resulting expressions for Se would 
decrease like exp(— A) and exp(— 3A), respectively, where A is proportional to (r p /a) 3//2 3> 1. Therefore we shall ignore at 
the outset the terms in exp(±3inf), which simply means that the method of averaging applies to them with greater accuracy 
because of their faster frequency. At this stage, then, we have 



. _ m :i a 3 n 
6 ~ Mi^R 5 " 



( a(iR 2 + -it 2 + -(ft - -7T7 2 ) + e mt h{R - -C + + ^ ) 



(A8) 



here, as before, we have written R = £a + 77b + £c, where c = a x b. 
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h 



When we integrate with respect to t we require integrals of the form 



■ exp(int) dt, 



(A9) 



where X, and Xj are two of the non-zero components of R = XA + Y~B. Using two standard equations for hyperbolic motion, 
viz. 



n't — e sinh F — F, 
and 

R — die coshF — 1), 
where 

n a — GM123, 
we write this as 



h = 



• exp 



— (e sinh F — F) 

1' 



dF. 



(A10) 
(All) 
(A12) 

(A13) 



Now we recall the assumption that the encounter is slow (eq. [2]), which implies that the exponent here is rapidly 

oscillating compared with the time scale of variations of the other factors in the integrand. The asymptotic treatment of 

such integrals is standard. The idea is to evaluate the integral approximately by deforming the contour of integration off 

the real axis to the saddle point of the exponent ("method of steepest descents"), which occurs in this case at R = 0, by 

equation (All). (The process is rather like that adopted by Heggie 1975 in evaluating the change in energy, and the relevant 

contour is illustrated there.) The rest of the integrand in equation (A13) is singular here, however, and so we first reduce the 

dR 

singularity using partial integration. From equation (All) we have — — = a'e smhF, and so equation (A13) becomes 

dF 



T 1 [°° 1 d ( XiXj r. n , -n 

h = 5— / ttt — - < — — i- exp 1 — (e sinh F — F) \ } dF, 

3n'a' 2 e' R 3 dF I sinh F l [n' y '\ J 

o — OC 



(A14) 



though we have already had to deform the contour to avoid the singularity at F = 0. In performing the differentiation, to 
leading order we need only include the derivative of the exponential, since the factor n/n' brings down a factor of order 
(r-p/a) 3 / 2 . To leading order we therefore have 



h = 



3GM 123 e 



J — t 



XiXj 
R 2 sinhF 



exp 



-j (e sinh F — F) 



dF, 



where we have made use of equations (All) and (A12). 
A second application of the same trick leads to 



i 2 r 

ja' 2 n'e' 2 J_ ( 



XiXi 



■ exp 



i^-(e smhF - F) 
n' 



3GM 12 3a' 2 n'e' 2 J sinh 2 F 

and now we turn attention to the saddle point at R = 0, i.e 



dF, 



F — i arccos(l/e'). In its vicinity we have 



e sinh F — F ~ i\J e' 2 — 1 — i arccos — + ^ (F — i arccos — ) 2 \/ e' 2 — 1. 

e' 2 e 

Now the integration is trivial, since in the rest of the integrand we can substitute values at the saddle, and so 



2tt m; 



■ exp 



/ M12 
VM123 



ray 



' 2 \/e' 2 - 1 - arccos(l/e') 
(e' - I) 3 / 2 



3 G 1 /2M 1 5 2 / 3 4 a 9 / 4 rp /4 (e' + 1)5/4 
though we have yet to substitute values for Xi and Xj at the saddle point, where 
r p (e' + 1) 



R = 



-(A + tB). 

c 

Now we return to equation (A8). Using equations (A18) and (A19) we readily find that 

V2 /r „ X 3/2 ^/ e ,2 _ 1 _ arccos (l/ e ') 



(A15) 



(A16) 



(A17) 



(A18) 



(A19) 



W + iy 



3M 1 5 2 / 30 3 / 4 e' : 



x exp 



/ Mia \ 
V M123 / 



3» 9 ^ ^ ~ ~ 15-~»- 

-a-Aa B+-(a Ab-A-a Bb B) + —h Ab B 



(e' - l) 3 /2 
+ 



(A20) 



+ b 



-^(a.A) 2 + ^(a.B) 2 -^(a.Ab.B + a.Bb.A) + ^(b.A) 2 -^b.B) 2 
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The scalar products here are easily evaluated in terms of the usual Keplerian elements u, 0, and i, where fi is measured 
along the plane of motion of the binary from the direction of a to the ascending node, and u> is measured from the node to 
the direction of A. The last (geometric) factor in expression (A20) reduces, after some lengthy expressions, to the surprisingly 
compact expression 



3 9 
— -sin2tjsin i + — sin 2(fi + ui) (1 + cos i) 



+ b 



3 9 
--cos2cjsin i — — cos 2(0, + uj) (1 + cos i) 



(A21) 



whence the magnitude of the vector e is given by equation (12). When e' = 1 (parabolic motion) the exponential factor is 

— \ C^f 1 ' 2 ~) ( — V , and so we have arrived at equation (13). 
3 V M123 J \ a J 

A4 Initially Circular Binaries: The Impulsive Limit 

We consider briefly the case in which the inequality in equation (2) is reversed. The duration of the encounter is then very 
short compared with the binary period. This case appears unimportant for most astrophysical applications, and we include it 
here mainly for the discussion in Section 3.2. 

In equation (6) we model F as a delta function, and then this equation integrates to 

GMi 2 e = 2(<5r.r)r - (<5r.r)r, (A22) 

where 8i is the change in relative velocity of the components due to the encounter. (The second term in eq. [6] vanishes for 
a circular binary.) We evaluate <5r by integrating the quadrupole term of equation (4), keeping r fixed and taking for R the 
form corresponding to a rectilinear orbit. In the notation of Appendix Al this is R = r p A + VfB, if t = at the time of 
pericentre. It follows that 

Sh = (4AiAj + 2B l B J - 25a), (A23) 

where summation over j is implied, and the 5 on the right is the Kronecker symbol. 

When this result is substituted into equation (A22) we also write r = oa and r = nab, where the unit vectors have the 
meaning assigned to them in §A2, and n is defined in equation (Al). It quickly follows that 

GM 12 e = G ™* n f { 4 (2a AbA + aBbB)a-2 [2 (a ■ A) 2 + (a ■ B) 2 - l] £>} . (A24) 
Vr p 

The two terms are orthogonal, and so e is easily found. 

In the case of coplanar motion, we introduce the angle 0, defined to be the angle between the vectors a and A. Then the 
direction cosines are easily expressed in terms of 9, and equation (28) follows easily. 



APPENDIX B: THE EQUAL-MASS CASE 

The purpose of this appendix is to show that, for a circular binary with components of equal mass m, the induced eccentricity 
vanishes to all orders of the perturbation parameter a/r p . In first order perturbation theory this is easy to understand. All 
terms in equation (4) with n odd contribute nothing if mi = mi. Hence F contains only odd monomials in the components 
of the vector r, and then it is easy to see that the right side of equation (6) vanishes after averaging over one binary period, 
if e = 0. The following argument extends this reasoning to perturbation theory at all orders. 

Higher-order perturbation theory is best addressed using Hamiltonian methods, and here we describe the appropriate 
variables in the simplest non-trivial case: coplanar motion of all three stars. The relative motion of the binary components is 
conveniently handled in action-angle variables, which, for planar Keplerian motion, can be taken as equivalent to the classical 
Delaunay variables L, G, I and g (cf. Plummer 1918). The variables L and G are expressible in terms of the semi-major axis 
a and eccentricity e by means of 

a = 2L 2 /(Gm 3 ) (Bl) 
and 

G = Ly/l-e 2 . (B2) 

The conjugate variable I is the mean anomaly, which increases uniformly with time in unperturbed motion, passes the value 
at pericentre, and increases by 2tt in each revolution. The variable g is the longitude of the pericentre, which we assume is 
measured from the direction of one axis of an inertial frame. The Hamiltonian function for unperturbed Keplerian motion is 
a function of L alone, denoted by Hb(L). 

The motion of the third star we continue to describe by its position vector R with respect to the barycentre of the 
binary, and we introduce the conjugate momentum P. Then the Keplerian approximation to its motion can be described by 
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a Hamiltonian function H3 (R, P). 

The whole Hamiltonian can now be written as 

H = H b (L) + H 3 (R, P) + 5ft(r(L, G, I, g), R), (B3) 

where the three terms on the right correspond, respectively, to the binary, the Keplerian motion of the third body, and the 
perturbation. The last of these terms is, except for the terms with n = and 1 and a factor involving the masses, nothing 
other than the expression whose gradient appears in equation (4), but we write r in terms of the Delaunay variables. 

In the case of equal masses, 5ft is invariant under the transformation r — > — r, i.e. the transformation g — > g + -k, with the 
other three canonical variables unchanged. To see that these transformations are equivalent, it is useful to consider the case 
of a binary of low eccentricity. Approximate expressions for the relative position vector r of the two components are easily 
found (cf. Plummer 1918 again) as 

3 1 1 

r = a(cos/ — -e + -ecos2/, sin / + -esin2/), (B4) 

if the components are resolved along the directions of the major and minor axes of the relative orbit, and 

3 1 3 1 

r = a(cos(g + I) — -ecos g + -ecos(2Z + g), sin(g + /) — -esin <; + -esin(2/ + g)), (B5) 

relative to the same frame as that in which g is measured. The effect of the transformation g — > g + tv is obvious. 

Because g is ill-defined when e becomes very small, it is better to use the Poincare variables L' , I' , £ and 77, defined by 
L' — L, 1' = l + g, 

(B6) 




r/ = -V2G'sing (B7) 
where 

G' = L-G~ (l/2)e 2 L, (B8) 

by equation (B2). These are still canonical, I' and 77 being conjugate to L' and £, respectively. In these variables equation (B5) 
becomes 

r = a( cos I' - 5 JL + _L= (£ cos 2 l' - 77 sin 2/') , sin l' + 5 -^= + — ^= (£ sin 21' + r) cos 2/') ) (B9) 
2^1/ 2VX 7 2 *s/L> 2VJ7 

approximately, and the transformation I' — > I' + tv, £ — > — £, 77 — » —77 sends r — » — r. 

In general, r is a 27r-periodic function of I' , and can be expanded as a Fourier series r = ^ n r„(L', £, 77) exp(mi'), where 
the coefficients r n are actually power series in £ and 77. This is illustrated at lowest non-trivial order by equation (B9). Similarly, 
we may expand the perturbing function as a series 5ft = X)^°oo a ™(L' ' , V> R) exp(inZ'), where again the coefficients are power 
series in the variables £ and 77. Now the invariance of the Hamiltonian under the transformation r — > — r (for equal masses) 
implies that this series is invariant under the transformation /'—»/' + 7r, £ — » — £, 77 — » —77. If the coefficient a„ contains a 
term ^rj 13 it follows that a + (3 is even if n is even, and a + (3 is odd if n is odd. 

The aim of Hamiltonian perturbation theory is to remove those terms which depend on quickly- varying angles. At lowest 
order, this is equivalent to averaging over the fast motions. Since I or /' is the one rapidly changing quantity, this leaves only 
those terms in the Fourier series which have n — 0. These terms have as coefficients power series in £ and 77 which have even 
"order" (i.e. the value of a + f3). When these are differentiated with respect to these variables (to give Hamilton's equations 
for their evolution) the resulting expression has terms of only odd order in these variables. It follows, therefore, that one 
solution is £ = 77 = 0, and since this solution yields the correct initial condition for a circular binary, it follows that the 
induced eccentricity vanishes (in first order perturbation theory). 

All we have to do now is to show that terms of the wrong parity will not be introduced in the perturbation procedure. In 
the standard procedure, which is described in many texts (e.g. Born 1960; Goldstein 1980; Lichtenberg & Liebermann 1983), 
the terms with n =fc are removed by a near-identity canonical transformation. For example, a term a n exp(ini') is removed 
by a transformation to new variables (denoted here by stars) in which 

I -I +Q U7 e , £ = £ +—e , (BIO) 

etc., where 

S 1 (77, R, L'\C, P*) = g'g ^i* ) ' P * ) (Bll) 
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(cf. Lichtenberg & Lieberman 1983, eq. [2.2.17]). 

Now it is clear from inspection of these equations (and the corresponding equations for the transformation of the other 
variables), along with our result on the parity of the coefficient a n , that the transformation I' — > I' + n, £ — > — £, r\ — > —77 
corresponds, in the new variables, to the transformation /'* — > V* + n, — > — £*, ?7* — » — r;*. Since the Hamiltonian was 
invariant under the transformation in the old variables, it will also be invariant under the transformation in the new ones. 
It follows that, when the Hamiltonian is expressed in terms of the new variables, the new coefficients a n will have the same 
parity property as the old ones. It follows that the induced eccentricity will again vanish for an initially circular binary. By 
induction, this conclusion follows for perturbation theory at all orders. 



